%% Figures

%% Figure 1: Changes in Cutoffs Over Time

head(Options)
programs=Options(:,{'proceso','demrecode','punt_lastAdmit','punt_last_waitAdmit','not_filled','slots','MajorName'});

listPrograms=table2array(unique(Options(:,{'demrecode'})));
listPrograms(isnan(listPrograms))=[];

listPrograms(listPrograms>=38000)=[];

Cutoffs=NaN(size(listPrograms,1),3);
NotFilled=Cutoffs;
Slots=Cutoffs;
WaitListCutoffs=Cutoffs;
for year=2010:2012

    for j=1:size(listPrograms,1)
        p=listPrograms(j);
        
         loc= find(programs.demrecode==p &     programs.proceso==year);
         if ~isempty(loc)
         Cutoffs(j,year-2009)=programs.punt_lastAdmit(loc(1));
         WaitListCutoffs(j,year-2009)=programs.punt_last_waitAdmit(loc(1));
         NotFilled(j,year-2009)=programs.not_filled(loc(1));
         Slots(j,year-2009)=programs.slots(loc(1));
         end
    end
end

dif12=Cutoffs(:,2)-Cutoffs(:,3);
pick=Slots(:,2)>1 & Slots(:,3)>1 & ~isnan(dif12);

sum(dif12(pick)>2500)/sum(pick)
sum(dif12(pick)>5000)/sum(pick)

pick=Slots(:,2)>1 & Slots(:,3)>1 & ~isnan(dif12);
sum(dif12(pick)>2500 & Cutoffs(pick,2)<60000)/sum(pick & Cutoffs(:,2)<60000)
sum(dif12(pick)>2500 & Cutoffs(pick,2)>60000)/sum(pick & Cutoffs(:,2)>60000)

sum(dif12(pick)>5000 & Cutoffs(pick,2)<60000)/sum(pick & Cutoffs(:,2)<60000)

dif10=Cutoffs(:,2)-Cutoffs(:,1);
pick=Slots(:,2)>1 & Slots(:,3)>1 & ~isnan(dif10);
sum(dif10(pick)>2500)/sum(pick)
sum(dif10(pick)<-2500)/sum(pick)

sum(dif10(pick)>2500 & Cutoffs(pick,2)<60000)/sum(pick & Cutoffs(:,2)<60000)
sum(dif10(pick)>2500 & Cutoffs(pick,2)>60000)/sum(pick & Cutoffs(:,2)>60000)


sum((Cutoffs(:,2)-Cutoffs(:,3))>0)/sum((pick))
sum((Cutoffs(:,1)-Cutoffs(:,2))>0)/sum((pick))




dif=abs(Cutoffs(:,2)-Cutoffs(:,3));
pick=Slots(:,2)>1 & Slots(:,3)>1;

sum(dif(pick)>2500)/sum((pick))
sum(dif(pick)>5000)/sum((pick))
sum(dif(pick)>7500)/sum((pick))

dif=abs(Cutoffs(:,2)-Cutoffs(:,1));
pick=Slots(:,2)>1 & Slots(:,3)>1;

sum(dif(pick)>2500)/sum((pick))
sum(dif(pick)>5000)/sum((pick))
sum(dif(pick)>7500)/sum((pick))


C=parula(10);
pick=NotFilled(:,2)==0 & NotFilled(:,3)==0 & Slots(:,2)>20 & Slots(:,3)>20;





sum(pick)
figure(1)
ax=subplot(1,2,2)


% Calculate centerline
Bands=[];
centerline = 45000:100:80000;
Bands(:,1)=centerline-7500;
Bands(:,2)=centerline-5000;
Bands(:,3)=centerline-2500;
Bands(:,4)=centerline+2500;
Bands(:,5)=centerline+5000;
Bands(:,6)=centerline+7500;
Bands = Bands/100;
centerline = centerline/100;

hold on
xplot = centerline([1:end end:-1:1])';

col=[C(6,:);  C(5,:); C(4,:); C(4,:);  C(5,:); C(6,:);]

for i = 1:3
    ph(i) = patch(xplot, [Bands(:,i);  flipud(Bands(:,end-i+1))]', col(i,:) ,'EdgeColor', col(i,:), 'FaceAlpha', 0.3);
end
lh = line(centerline, centerline, 'LineWidth', 2, 'Color', C(1,:));
ph = flipud(ph);
ax1=gca
xlim([450 800]);


h=scatter(Cutoffs(pick,2)/100,Cutoffs(pick,3)/100,10    ,'MarkerFaceColor',[0.75 0.75 0.75],'MarkerEdgeColor',C(2,:));
grid on

statsbajo = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,3)/100,:);
statsbajo25 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,3)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,3)/100 <= 25);
statsbajo50 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,3)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,3)/100 > 25 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,3)/100 <= 50 );
statsbajo75 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,3)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,3)/100 > 50 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,3)/100 <= 75 );

length(statsbajo) % 445 programas 2012<2011 
length(statsbajo25) % 335 programas 2012<2011 
length(statsbajo50) % 94 programas 2012<2011 
length(statsbajo75) % 13 programas 2012<2011 
        
        
        
        
statsbajo = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,1)/100,:);
statsbajo25 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,1)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,1)/100 <= 25);
statsbajo50 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,1)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,1)/100 > 25 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,1)/100 <= 50 );
statsbajo75 = Cutoffs(Cutoffs(pick,2)/100 > Cutoffs(pick,1)/100 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,1)/100 > 50 & ...
            Cutoffs(pick,2)/100-Cutoffs(pick,1)/100 <= 75 );

length(statsbajo) % 445 programas 2012<2011 
length(statsbajo25) % 335 programas 2012<2011 
length(statsbajo50) % 94 programas 2012<2011 
length(statsbajo75) % 13 programas 2012<2011 
        



ylim([450 800]);
xlim([450 800]);

grid on 
ax1=gca
xlabel('Cutoff Score in 2011 (G25)','Interpreter','latex')
ylabel('Cutoff Score in 2012 (G25)','Interpreter','latex')
ax1.GridAlpha=0.1;
box on
set(gca,'TickLabelInterpreter','latex')
title('Comparison 2012','Interpreter','latex')

loc= find(abs(Cutoffs(:,2)-Cutoffs(:,3))>7500 & pick==1);


pickName=(ismember(programs.demrecode,listPrograms(loc)) & programs.proceso==2012)
programs.MajorName(pickName)

% Cutoffs vs WaitlistCutoffs 
C=parula(10);
pick=NotFilled(:,1)==0 & NotFilled(:,2)==0 & Slots(:,2)>20 & Slots(:,3)>20;
sum(pick)
ax=subplot(1,2,1)



% Calculate centerline
Bands=[];
centerline = 45000:100:80000;
Bands(:,1)=centerline-7500;
Bands(:,2)=centerline-5000;
Bands(:,3)=centerline-2500;
Bands(:,4)=centerline+2500;
Bands(:,5)=centerline+5000;
Bands(:,6)=centerline+7500;
Bands = Bands/100;
centerline = centerline/100;

hold on
xplot = centerline([1:end end:-1:1])';

col=[C(6,:);  C(5,:); C(4,:); C(4,:);  C(5,:); C(6,:);]

for i = 1:3
    ph(i) = patch(xplot, [Bands(:,i);  flipud(Bands(:,end-i+1))]', col(i,:) ,'EdgeColor', col(i,:), 'FaceAlpha', 0.3);
end
lh = line(centerline, centerline, 'LineWidth', 2, 'Color', C(1,:));
ph = flipud(ph);
ax1=gca
xlim([450 800]);

h=scatter(Cutoffs(pick,2)/100,Cutoffs(pick,1)/100,10    ,'MarkerFaceColor',[0.75 0.75 0.75],'MarkerEdgeColor',C(2,:));
grid on


ylim([450 800]);
xlim([450 800]);

grid on 
ax1=gca
xlabel('Cutoff Score in 2011 (G25)','Interpreter','latex')
ylabel('Cutoff Score in 2010 (G25)','Interpreter','latex')
ax1.GridAlpha=0.1;
box on


loc= find(abs(Cutoffs(:,2)-Cutoffs(:,3))>7500 & pick==1);
set(gca,'TickLabelInterpreter','latex')
title('Comparison 2010','Interpreter','latex')
set(gcf, 'Units', 'Inches', 'Position', [1, 1, 13,  6], 'PaperUnits', 'Inches', 'PaperSize', [7.25, 9.125])


% Changes in Program Cutoffs Over Time
savefig(fullfile(pathOutput, 'figures',"cutoff_figure.fig")); 
saveas(gcf, fullfile(pathOutput, 'figures',"cutoff_figure.png")); 
exportgraphics(gcf, fullfile(pathOutput, 'figures','cutoff_figure.eps'),'BackgroundColor','none','ContentType','vector')

%% Figure 2: Diagram of the on-platform application process
% (HTML source here)
% <html>
% <body>
%  <script type="text/javascript" src="https://www.gstatic.com/charts/loader.js"></script>
% 
% <div id="sankey_multiple" style="width: 600px; height: 400px;"></div>
% 
% <script type="text/javascript">
%   google.charts.load("current", {packages:["sankey"]});
%   google.charts.setOnLoadCallback(drawChart);
%    function drawChart() {
%     var data = new google.visualization.DataTable();
%     data.addColumn('string', 'From');
%     data.addColumn('string', 'To');
%     data.addColumn('number', 'Students');
%     data.addRows([
%        [ 'Test-takers (250k)', 'Eligible (170k)', 169332 ],
%        [ 'Test-takers (250k)', 'Ineligible (80k)', 81426 ],
%        [ 'Eligible (170k)', 'Applied (85k)', 84512 ],
%        [ 'Eligible (170k)', 'Did not apply (85k)', 84820 ],
%        [ 'Ineligible (80k)', 'Applied (85k)', 1 ],
%        [ 'Applied (85k)', 'G25 Admission (68k)', 67803 ],
%        [ 'Applied (85k)', 'G25 Waitlist (17k)', 16709 ],
%        [ 'Did not apply (85k)', 'G25 Admission (68k)', 1 ],
%        [ 'G25 Admission (68k)', 'G25 Regular Enrollment (46k)', 45967 ],
%        [ 'G25 Admission (68k)', 'G25 Waitlist Enrollment (4k)', 3920 ],
%        [ 'G25 Admission (68k)', 'G8 Enrollment (8k)', 4352 ],
%        [ 'G25 Waitlist (17k)', 'G25 Waitlist Enrollment (4k)', 1679 ],
%        [ 'G25 Waitlist (17k)', 'G8 Enrollment (8k)', 3555 ]
%     ]);
% 
%     // Parula colors
%     var colors = ['#3E27A9', '#4847EC', '#3F70FF', '#2897EC',
%               '#08B5D1', '#32C6A0', '#81CC59', '#DDBD29',
%               '#FDD031', '#FAFB15', '#885826'];
% 
%     // Set chart options
%     var options = {
%       width: 800,
%       sankey: {
%         link: { color: { fill: '#D3D3D3', fillOpacity: 0.1 } },
%         node: { colors: colors,
%           label: { fontName: 'Times-Roman',
%                            fontSize: 14 } },
%         }
% 
%     };
% 
%     // Instantiate and draw our chart, passing in some options.
%     var chart = new google.visualization.Sankey(document.getElementById('sankey_multiple'));
%     chart.draw(data, options);
%    }
% </script>
% </body>
% </html>



%% Figure 3: Total Slots, Excess Slots, Program Yield
Students.Properties.VariableNames'
Students(1:10,:)
tabulate(Students.proceso)
tabulate(Students.current_cohort(Students.proceso==2011))
tabulate(Students.current_cohort(Students.proceso==2012))
tabulate(Students.current_cohort(Students.proceso==2010))

tabulate(Students.private_hs)
tabulate(Students.male)

  
TableA=[];TableB=[];TableC=[]; scoregrid=[];
Years=[2010:2012];inputGraphs={}; it=1;indexGraphs=[];
for y=1:3
    
    year_j=Years(y);
    pickYear=Apps.proceso==year_j;
    
    Enroll_Year=Enrollment(Enrollment.proceso==year_j,:);
    
    ListCodes=unique(Apps.cod(pickYear));
    for j=1:length(ListCodes)
        
        cod_j=ListCodes(j);

        enroll_j=Enrollment(Enrollment.proceso==year_j & cod_j==Enrollment.cod,:);
        
        apps_j=Apps(Apps.proceso==year_j & cod_j==Apps.cod,:);
        
        apps_j=sortrows(apps_j,{'lugar'});
        mat_j=ismember(apps_j.rut,enroll_j.rut);
        
        
        students=Students(Students.proceso==year_j,{'rut','bea','private_hs','male'});
        apps_j=join(apps_j,students,'Keys',{'rut'});
        
        
        
        waitlistMat=(apps_j.sit==25 & mat_j==1 & apps_j.punt_lastAdmit-apps_j.punt<5000 & apps_j.bea==0);
        waitlistMatBEA=(apps_j.sit==25 & mat_j==1 & apps_j.punt_lastAdmit-apps_j.punt<5000 & apps_j.bea==1);
        
        sit_j=[sum(apps_j.sit==26) sum(apps_j.sit==24) sum(apps_j.sit==25) sum(apps_j.sit==24 & mat_j==1) sum(waitlistMat) sum(waitlistMatBEA)];
        
        yeild_j=sum(apps_j.sit==24 & mat_j==1)/sum(apps_j.sit==24);
        
        classSize_j=[sum(mat_j==1)/sum(apps_j.sit==24) sum(mat_j==1)/apps_j.vacantes(1)];
        
        
        backdoor_j=1-sum(mat_j)/size(enroll_j,1);
        
        
        waitlist_j=sum(waitlistMat)/sum(mat_j);
        
        if waitlist_j>0
            
            if max(apps_j.lugar(waitlistMat))==min(apps_j.lugar(waitlistMat))
                waitlistCall_j=1;      
            else
                waitlistCall_j=sum(waitlistMat)/(max(apps_j.lugar(waitlistMat))-min(apps_j.lugar(waitlistMat)));
            end
            
            
            LastAdmit=(max(apps_j.lugar(waitlistMat)));
            cutoff=max(apps_j.lugar(apps_j.sit==24));
            
            xgrid=[[sort(cutoff:-10:5)-5]'; sort(cutoff+5:10:LastAdmit)']; 
            ygrid=[];
            scoregrid=[];
            jj=1;
            for jj=1:length(xgrid)
                ii=xgrid(jj);
                pickAbove=(apps_j.lugar<=ii+5 & apps_j.lugar>ii-5 & apps_j.sit==24 & ii<=cutoff );
                pickBelow=(apps_j.lugar<=ii+5 & apps_j.lugar>ii-5 & apps_j.sit==25  & ii>cutoff & apps_j.bea==0);
                
                ngrid(jj,1)=sum(pickAbove);
                ngrid(jj,2)=sum(pickBelow);
                
                ygrid(jj,1)=mean(mat_j(pickAbove==1 ));
                ygrid(jj,2)=mean(mat_j(pickBelow==1));
                
                scoregrid(jj,1)=mean(apps_j.punt(pickAbove==1 | pickBelow==1));
            end
            
            indexGraphs(it,1:2)=[year_j cod_j];
            inputGraphs{it,1}=[xgrid ygrid scoregrid];
            it=it+1;
        else
            waitlistCall_j=NaN;
        end
    
        TableA=[TableA ; year_j cod_j yeild_j waitlist_j waitlistCall_j backdoor_j classSize_j ];
        TableB=[TableB ; year_j cod_j sum(apps_j.sit==24 | (apps_j.sit==25 & apps_j.bea==0))>apps_j.vacantes(1) sum(apps_j.sit==24 | (apps_j.sit==25 & apps_j.bea==0))>apps_j.posted_slots(1)  sum(apps_j.sit==24 | mat_j==1)>apps_j.vacantes(1)   sum((apps_j.sit==25 & apps_j.bea==0) & mat_j==1)>0   ];
        TableC=[TableC ; year_j cod_j sit_j  ];
        
        
    end
end

TA = array2table(TableA,'VariableNames',{'Year','cod','AcceptedOffers','PercentMatFromWaitlist','ProbMatFromWaitlist','PercentMatFromOther','RatioMatOverNumberOffers','RatioMatOverNumberDesired'});
% TA.Properties.VariableNames'
TB = array2table(TableB,'VariableNames',{'Year','cod','ExcessApps2DesiredSlots','ExcessApps2PostedSlots','OffersAcceptedOverDesiredSlots','Went2Waitlist'});
% TB.Properties.VariableNames'
TC = array2table(TableC,'VariableNames',{'Year','cod','AppsAcceptedAtHigherRanked','AppsOffered','AppsWaitlisted','OffersAccepted','WaitlistOffersAccepted','WaitlistBEAOffersAccepted'});
% TC.Properties.VariableNames'




Options.Properties.VariableNames'

nansum(Options.seats(Options.proceso==2011))
nansum(Options.seats(Options.proceso==2012))
nansum(Options.slots(Options.proceso==2012))

Options(1:10,:)
sum(Options.slots(Options.proceso==2012))-sum(Options.posted_slots(Options.proceso==2011))


sum(Options.slots(Options.proceso==2011 & Options.Univ>=38))


% Offers made 
sum(Options.posted_slots(Options.proceso==2011))

sum(TC.WaitlistOffersAccepted)

pickYears=(TA.Year<2012);

% options with not enough applications to fill slots desired
sum(TB.ExcessApps2DesiredSlots(TB.Year==2011)==0)/sum(TB.Year==2011)
notEnough=(TB.ExcessApps2DesiredSlots==0)

% share of options that had enough apps to go to the waitlist 
sum(TB.Went2Waitlist(TB.Year==2011 & notEnough==0)==1)/sum(TB.Year==2011 & notEnough==0);

sum(TB.ExcessApps2DesiredSlots(TB.Year==2011)==1)

sum(TB.Went2Waitlist)


sum(TC.WaitlistBEAOffersAccepted(TC.Year==2011))

sum(TC.WaitlistOffersAccepted(TC.Year==2011))/(sum(TC.OffersAccepted(TC.Year==2011))+sum(TC.WaitlistOffersAccepted(TC.Year==2011)));


sum(TC.WaitlistOffersAccepted(TC.Year<2012))/(sum(TC.OffersAccepted(TC.Year<2012))+sum(TC.WaitlistOffersAccepted(TC.Year<2012)));


sum(TC.AppsOffered(pickYears(TC.Year==2011)))


TA.AcceptedOffers

sum(TA.PercentMatFromWaitlist>0)


sum(TC.WaitlistOffersAccepted(pickYears)>0)/sum(pickYears)

sum(TC.WaitlistOffersAccepted(pickYears))/sum(pickYears)


sum(Options.sobrecupo(Options.proceso==2011 & Options.Univ<38))
sum(Options.slots(Options.proceso==2011 & Options.Univ<38))

xbins=[0:.1:1];


figure(3)

cupos=Options.slots(Options.proceso==2011 & Options.Univ<38);
sobrecupos=(Options.sobrecupo(Options.proceso==2011 & Options.Univ<38))

prctile(sobrecupos,95)
prctile(cupos,95)
prctile(cupos,5)


subplot(3,2,1)
xbins=linspace(5,305,11);
[xi,bins]=hist(cupos,xbins);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(1,:),'EdgeColor','none')
xticks(xbins)
xlim([0 xbins(end)+10])
ylim([0 0.4])
yticks(linspace(0,.4,11))
xlabel('Total Slots ')
ylabel('Share of Programs')
set(gca,'TickLabelInterpreter','latex')
title('Distribution of Total Slots at Programs')
grid on



subplot(3,2,2)
xbins=linspace(0,100,11);
[xi,bins]=hist(sobrecupos,xbins);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(3,:),'EdgeColor','none')
xticks(xbins)
xlim([0 105])
ylim([0 0.4])
yticks(linspace(0,.4,11))
xlabel('Extra Slots')
ylabel('Share of Programs')
set(gca,'TickLabelInterpreter','latex')
title('Distribution of Extra Slots')
grid on


subplot(3,2,3)
xbins=0:.1:1;
ratio=Options.sobrecupo(Options.proceso==2011 & Options.Univ<38)./Options.slots(Options.proceso==2011 & Options.Univ<38);
prctile(ratio,95)
[xi,bins]=hist(ratio,xbins);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(5,:),'EdgeColor','none')
xticks(xbins)
xlim([0 1.05])
ylim([0 0.3])
yticks(linspace(0,.3,11))
xlabel('Ratio Extra Slots / Desired Slots ')
ylabel('Share of Programs')
set(gca,'TickLabelInterpreter','latex')
title('Distribution of Ratio Extra Slots  / Desired Slots')
grid on



subplot(3,2,4)
pickYears=(TA.Year==2011);
xbins=0:.1:1;
[xi,bins]=hist(TA.AcceptedOffers(pickYears),11);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(6,:),'EdgeColor','none')
xlim([0 1])
grid on
ylim([0 0.25])
ylabel('Share of Programs')
xlabel('Share of Initial Offers Acceptd (Yield)')
set(gca,'TickLabelInterpreter','latex')
title('Yield on Offers')
xticks(xbins)
xlim([0 1])
ylim([0 0.3])
yticks(linspace(0,.3,11))
grid on



subplot(3,2,5)
xbins=[0:.2:2];
[xi,bins]=hist(TA.RatioMatOverNumberDesired(pickYears),11);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(7,:),'EdgeColor','none')
xlim([0 2])
ylim([0 0.4])
xlabel('Final Matriculation / Desired Slots ')
ylabel('Share of Programs')
set(gca,'TickLabelInterpreter','latex')
title('Matriculation / Desired Slots ')
xticks(xbins)
xlim([0 2])
ylim([0 0.4])
yticks(linspace(0,.4,11))
grid on

subplot(3,2,6)
xbins=[0:.05:.5];
waitOfferShare=(TC.WaitlistOffersAccepted)./((TC.OffersAccepted)+(TC.WaitlistOffersAccepted));
ratio=TA.PercentMatFromWaitlist
[xi,bins]=hist(ratio(pickYears==1),xbins);

sum(pickYears==1 & ratio>0)
sum(pickYears==1 & ratio==0)/sum(pickYears==1)

[xi,bins]=hist(waitOfferShare(pickYears==1 & ratio>0),xbins);
bar(bins,xi/sum(xi),'LineWidth',1,'FaceColor',ParulaColors(8,:),'EdgeColor','none')
ylabel('Share of Programs')
xlabel('Share of Matriculation from Waitlist')
set(gca,'TickLabelInterpreter','latex')
title('Waitlist Matriculation / Total Matriculation')
xticks(xbins)
xlim([0 0.55])
xlim([0 0.5])
yticks(linspace(0,.5,11))
grid on

set(gcf, 'Units', 'Inches', 'Position', [1, 1, 9, 10], 'PaperUnits', 'Inches', 'PaperSize', [7.25, 9.125])

% Total Slots, Excess Slots, Program Yield
savefig(fullfile(pathOutput, 'figures',"summarySlots.fig")); 
saveas(gcf, fullfile(pathOutput, 'figures',"summarySlots.png")); 
exportgraphics(gcf, fullfile(pathOutput, 'figures','summarySlots.eps'),'BackgroundColor','none','ContentType','vector')

%% Figure 4: Case Study: Enrollment Probability at Economics - University of Chile - 2010/2011


% Ingeneria Comercial - Uchile : 1142
Options(strcmp(Options.MajorName,'Ingeniería Comercial')==1 & Options.Univ==11,:)
case1=find((strcmp(Options.MajorName,'Ingeniería Comercial')==1 & Options.Univ==11 & Options.proceso==2011));

interestList=[Options.cod(case1)];
insterestNames={'Business and Economics - University of Chile'};

for i=1:length(interestList)
    
pickCOD=interestList(i); 
pickRows=(TA.Year<2012 & TA.cod==pickCOD);
mean(TA.AcceptedOffers(pickRows))
mean(TA.PercentMatFromWaitlist(pickRows))
TA.RatioMatOverNumberDesired(pickRows)
TableC(TableC(:,1)<2012 & TableC(:,2)==pickCOD,:)
Options(Options.cod==pickCOD,{'seats','sobrecupo'})

end


close all
cod_j=interestList;
    
    pick1=find((indexGraphs(:,1)==2010 & indexGraphs(:,2)==cod_j));    
    pick2=find((indexGraphs(:,1)==2011 & indexGraphs(:,2)==cod_j));

    Temp1=inputGraphs{pick1,1};
    Temp2=inputGraphs{pick2,1};
    
    l=min([size(Temp1,1) size(Temp2,1)]);
    
      xI=Temp1(1:l,1);
      yI=(Temp1(1:l,2:3)+Temp2(1:l,2:3))/2;
      sI=(Temp1(1:l,4)+Temp2(1:l,4))/2;
    
      figure(4)
      
      %yyaxis left
      bar(xI-4,yI(:,1),'LineWidth',1,'FaceColor',c(105,:),'EdgeColor','none')
      hold on
      bar(xI-4,yI(:,2),'LineWidth',1,'FaceColor',c(10,:),'EdgeColor','none')
      ylim([0 1])
      ylabel('Pr(Mat|Rank)')
      xlim([xI(1)-4 xI(end)])
      xlabel('Applicant Rank')
      grid on
set(gca,'TickLabelInterpreter','latex')   
set(gcf, 'Units', 'Inches', 'Position', [1, 1, 9, 3], 'PaperUnits', 'Inches', 'PaperSize', [7.25, 9.125])

% Case Study: Enrollment Probability at Economics - University of Chile - 2010/2011
savefig(fullfile(pathOutput, 'figures',"CaseStudyUChile.fig")); 
saveas(gcf, fullfile(pathOutput, 'figures',"CaseStudyUChile.png")); 
exportgraphics(gcf, fullfile(pathOutput, 'figures','CaseStudyUChile.eps'),'BackgroundColor','none','ContentType','vector')


%% Figure 5: Enrollment probabilities for G25 admits

G33adm=Apps(Apps.sit==24,{'rut','proceso','cod','punt'});
G33adm.Univ=floor(G33adm.cod./(10.^(-1+floor(log10(G33adm.cod)))));
G33admenroll=outerjoin(G33adm,unique(Enrollment(:,{'rut','proceso','Univ','dropout','grad6'})),'Type','left','Keys',{'rut','proceso'},'MergeKeys',true);
G33admenroll.G25enroll=(G33admenroll.Univ_G33adm==G33admenroll.Univ_right);
G33admenroll.G8enroll=(G33admenroll.Univ_G33adm<G33admenroll.Univ_right).*(G33admenroll.Univ_right>37);%.*(G33admenroll.proceso<2012);
G33admenroll.otherornot=isnan(G33admenroll.Univ_right);
G33admenroll.specialcase=1-G33admenroll.G25enroll-G33admenroll.G8enroll-G33admenroll.otherornot;
sumG33admenroll=grpstats(G33admenroll(:,{'proceso','G25enroll','G8enroll','otherornot','specialcase'}),'proceso')
G33admG25enr=sumG33admenroll.GroupCount.*sumG33admenroll.mean_G25enroll;
G33admG8enr=sumG33admenroll.GroupCount.*sumG33admenroll.mean_G8enroll;
G33admothernot=sumG33admenroll.GroupCount.*sumG33admenroll.mean_otherornot;
G33admspecial=sumG33admenroll.GroupCount.*sumG33admenroll.mean_specialcase;

G25admenroll=G33admenroll(G33admenroll.Univ_G33adm<38,:);
sumG25admenroll=grpstats(G25admenroll(:,{'proceso','G25enroll','G8enroll','otherornot','specialcase'}),'proceso')


figure(6)
b = bar(sumG25admenroll.proceso,[sumG25admenroll.mean_G25enroll, sumG25admenroll.mean_G8enroll,sumG25admenroll.mean_specialcase, sumG25admenroll.mean_otherornot],'stacked');
th = title("G25 Admits and Enrollment Probabilities");
set(b(1),'FaceColor',c(80,:))
set(b(2),'FaceColor',c(200,:))
set(b(3),'FaceColor',c(135,:))
set(b(4),'FaceColor', uint8([200 200 200]))
grid
ylabel('Percentage of students admitted to G25')
ylim([0.5 1])
xlim([2009.2 2012.8])
text(2009.65,0.6,'Assigned and') 
text(2009.65,0.58,'Enrolled G25') 
text(2010.65,0.6,'Assigned and') 
text(2010.65,0.58,'Enrolled G25') 
text(2011.65,0.6,'Assigned and') 
text(2011.65,0.58,'Enrolled G25') 
text(2009.7,0.75,'Enrolled G8') 
text(2010.7,0.75,'Enrolled G8')

text(2009.77,0.79,'Other G25') 
text(2010.77,0.787,'Other G25') 
text(2011.77,0.792,'Other G25') 

text(2009.85,0.9,'Other') 
text(2010.85,0.9,'Other') 
text(2011.85,0.9,'Other') 

yticks([0.5:0.1:1])
ytickformat('%.1f')

hold off
set(gca,'TickLabelInterpreter','latex') 
set(gcf, 'Units', 'Inches', 'Position', [1, 1, 8, 5], 'PaperUnits', 'Inches', 'PaperSize', [7.25, 9.125])
% Enrollment probabilities for G25 admits
saveas(gcf,fullfile(pathOutput, 'figures', 'enrollment_outcomes_G25adm.png'))

exportgraphics(gcf, fullfile(pathOutput, 'figures','enrollment_outcomes_G25adm.eps'),'BackgroundColor','none','ContentType','vector')

%% Figure 6: Enrollment and Dropout probability, G25
maxpsu=820;
minpsu=500;
G25admenrollscore=G25admenroll([G25admenroll.punt<=maxpsu*100].*[G25admenroll.punt>=minpsu*100]==1,{'proceso','punt','G25enroll'});
G25enrolldropgradscore=G25admenroll([G25admenroll.punt<=maxpsu*100].*[G25admenroll.punt>=minpsu*100].*G25admenroll.G25enroll==1,{'proceso','punt','dropout','grad6'});

binsize=70;
G25admenrollscore.score=floor(G25admenrollscore.punt/100);
G25admenrollscore.psu=G25admenrollscore.punt/100;
prenrollG25score=zeros(maxpsu-minpsu+1-binsize,height(sumG25admenroll));
G25admenrollscore2010=G25admenrollscore(G25admenrollscore.proceso==2010,{'score','G25enroll'});
G25admenrollscore2011=G25admenrollscore(G25admenrollscore.proceso==2011,{'score','G25enroll'});
G25admenrollscore2012=G25admenrollscore(G25admenrollscore.proceso==2012,{'score','G25enroll'});
for s=[minpsu:maxpsu]
    if s>=maxpsu-binsize+1
        break
    end
    prenrollG25score(s-minpsu+1,1)=mean(table2array(G25admenrollscore2012([G25admenrollscore2012.score>=s].*[G25admenrollscore2012.score<s+binsize]==1,'G25enroll')));
    prenrollG25score(s-minpsu+1,2)=mean(table2array(G25admenrollscore2011([G25admenrollscore2011.score>=s].*[G25admenrollscore2011.score<s+binsize]==1,'G25enroll')));
    prenrollG25score(s-minpsu+1,3)=mean(table2array(G25admenrollscore2010([G25admenrollscore2010.score>=s].*[G25admenrollscore2010.score<s+binsize]==1,'G25enroll')));
end


x = [minpsu:maxpsu-binsize];
xx = 500:20:750;
y1 = prenrollG25score(:,1);
y2 = prenrollG25score(:,2);
y3 = prenrollG25score(:,3);
yy1 = spline(x,y1,xx);
yy2 = spline(x,y2,xx);
yy3 = spline(x,y3,xx);
figure(71)
    hold on
    plot(xx,yy1,'Color', c(80,:),'LineWidth',2); 
    plot(xx,yy2,'Color', c(160,:),'LineWidth',2); 
    plot(xx,yy3,'Color', c(240,:),'LineWidth',2); 
    xlabel("Score-Bin (Floor, Range="+ num2str(binsize)+ ")")
    legend({'2012','2011','2010'},'Location','southeast','Interpreter','latex') 
    legend('boxoff')
    set(gca,'FontSize',12)
    grid on
    ytickformat('%.2f')
    box on;
    hold off
    % Enrollment probability, G25
    set(gca,'TickLabelInterpreter','latex') 
    saveas(gcf, fullfile(pathOutput, 'figures', 'enrollment_score.png'));
    exportgraphics(gcf, fullfile(pathOutput, 'figures','enrollment_score.eps'),'BackgroundColor','none','ContentType','vector')

close


%% Figure 7: Enrollment and dropout probability, conditional on scores, gender, and type of school

maxpsu=770;
minpsu=500;
G25admenroll=G33admenroll(G33admenroll.Univ_G33adm<38,:);
G25admenroll=outerjoin(G25admenroll,unique(Students(:,{'rut','male','private_hs'})),'Keys',{'rut'},'MergeKeys',true);
G25admenroll=rmmissing(G25admenroll,'DataVariables',{'proceso'});
G25admenrollscore=G25admenroll([G25admenroll.punt<=maxpsu*100].*[G25admenroll.punt>=minpsu*100]==1,{'proceso','punt','G25enroll','male','private_hs'});
binsize=70;
G25admenrollscore.score=floor(G25admenrollscore.punt/100);
G25admenrollscore.psu=G25admenrollscore.punt/100;
prenrollG25score1=zeros(maxpsu-minpsu+1-binsize,height(sumG25admenroll));
G25admenrollscore20101=G25admenrollscore([G25admenrollscore.proceso==2010].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
G25admenrollscore20111=G25admenrollscore([G25admenrollscore.proceso==2011].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
G25admenrollscore20121=G25admenrollscore([G25admenrollscore.proceso==2012].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
prenrollG25score2=zeros(maxpsu-minpsu+1-binsize,height(sumG25admenroll));
G25admenrollscore20102=G25admenrollscore([G25admenrollscore.proceso==2010].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
G25admenrollscore20112=G25admenrollscore([G25admenrollscore.proceso==2011].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
G25admenrollscore20122=G25admenrollscore([G25admenrollscore.proceso==2012].*[G25admenrollscore.male==1].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
prenrollG25score3=zeros(maxpsu-minpsu+1-binsize,height(sumG25admenroll));
G25admenrollscore20103=G25admenrollscore([G25admenrollscore.proceso==2010].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
G25admenrollscore20113=G25admenrollscore([G25admenrollscore.proceso==2011].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
G25admenrollscore20123=G25admenrollscore([G25admenrollscore.proceso==2012].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==1]==1,{'score','G25enroll'});
prenrollG25score4=zeros(maxpsu-minpsu+1-binsize,height(sumG25admenroll));
G25admenrollscore20104=G25admenrollscore([G25admenrollscore.proceso==2010].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
G25admenrollscore20114=G25admenrollscore([G25admenrollscore.proceso==2011].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
G25admenrollscore20124=G25admenrollscore([G25admenrollscore.proceso==2012].*[G25admenrollscore.male==0].*[G25admenrollscore.private_hs==0]==1,{'score','G25enroll'});
for s=[minpsu:maxpsu]
    if s>=maxpsu-binsize+1
        break
    end
    prenrollG25score1(s-minpsu+1,1)=mean(table2array(G25admenrollscore20121([G25admenrollscore20121.score>=s].*[G25admenrollscore20121.score<s+binsize]==1,'G25enroll')));
    prenrollG25score1(s-minpsu+1,2)=mean(table2array(G25admenrollscore20111([G25admenrollscore20111.score>=s].*[G25admenrollscore20111.score<s+binsize]==1,'G25enroll')));
    prenrollG25score1(s-minpsu+1,3)=mean(table2array(G25admenrollscore20101([G25admenrollscore20101.score>=s].*[G25admenrollscore20101.score<s+binsize]==1,'G25enroll')));
    prenrollG25score2(s-minpsu+1,1)=mean(table2array(G25admenrollscore20122([G25admenrollscore20122.score>=s].*[G25admenrollscore20122.score<s+binsize]==1,'G25enroll')));
    prenrollG25score2(s-minpsu+1,2)=mean(table2array(G25admenrollscore20112([G25admenrollscore20112.score>=s].*[G25admenrollscore20112.score<s+binsize]==1,'G25enroll')));
    prenrollG25score2(s-minpsu+1,3)=mean(table2array(G25admenrollscore20102([G25admenrollscore20102.score>=s].*[G25admenrollscore20102.score<s+binsize]==1,'G25enroll')));
    prenrollG25score3(s-minpsu+1,1)=mean(table2array(G25admenrollscore20123([G25admenrollscore20123.score>=s].*[G25admenrollscore20123.score<s+binsize]==1,'G25enroll')));
    prenrollG25score3(s-minpsu+1,2)=mean(table2array(G25admenrollscore20113([G25admenrollscore20113.score>=s].*[G25admenrollscore20113.score<s+binsize]==1,'G25enroll')));
    prenrollG25score3(s-minpsu+1,3)=mean(table2array(G25admenrollscore20103([G25admenrollscore20103.score>=s].*[G25admenrollscore20103.score<s+binsize]==1,'G25enroll')));
    prenrollG25score4(s-minpsu+1,1)=mean(table2array(G25admenrollscore20124([G25admenrollscore20124.score>=s].*[G25admenrollscore20124.score<s+binsize]==1,'G25enroll')));
    prenrollG25score4(s-minpsu+1,2)=mean(table2array(G25admenrollscore20114([G25admenrollscore20114.score>=s].*[G25admenrollscore20114.score<s+binsize]==1,'G25enroll')));
    prenrollG25score4(s-minpsu+1,3)=mean(table2array(G25admenrollscore20104([G25admenrollscore20104.score>=s].*[G25admenrollscore20104.score<s+binsize]==1,'G25enroll')));
end

x = [minpsu:maxpsu-binsize];
xx = 500:20:700;
y11 = prenrollG25score1(:,1);
y12 = prenrollG25score1(:,2);
y13 = prenrollG25score1(:,3);
y21 = prenrollG25score2(:,1);
y22 = prenrollG25score2(:,2);
y23 = prenrollG25score2(:,3);
y31 = prenrollG25score3(:,1);
y32 = prenrollG25score3(:,2);
y33 = prenrollG25score3(:,3);
y41 = prenrollG25score4(:,1);
y42 = prenrollG25score4(:,2);
y43 = prenrollG25score4(:,3);

yy11 = spline(x,y11,xx);
yy12 = spline(x,y12,xx);
yy13 = spline(x,y13,xx);

yy21 = spline(x,y21,xx);
yy22 = spline(x,y22,xx);
yy23 = spline(x,y23,xx);

yy31 = spline(x,y31,xx);
yy32 = spline(x,y32,xx);
yy33 = spline(x,y33,xx);

yy41 = spline(x,y41,xx);
yy42 = spline(x,y42,xx);
yy43 = spline(x,y43,xx);

figure(8)
    subplot(2,2,1)
    hold on
    set(gca,'TickLabelInterpreter','latex') 
    title('Male, Private')
    plot(xx,yy11,'Color', c(80,:),'LineWidth',1.2); 
    plot(xx,yy12,'Color', c(160,:),'LineWidth',1.2); 
    plot(xx,yy13,'Color', c(240,:),'LineWidth',1.2); 
    ylabel('Pr(enrollment)')
    ylim([0.5 1])
    grid on
    ytickformat('%.2f')
    box on;
    hold off
    subplot(2,2,2)
    hold on
    set(gca,'TickLabelInterpreter','latex') 
    title('Male, Public')
    plot(xx,yy21,'Color', c(80,:),'LineWidth',1.2); 
    plot(xx,yy22,'Color', c(160,:),'LineWidth',1.2); 
    plot(xx,yy23,'Color', c(240,:),'LineWidth',1.2); 
    ylabel('Pr(enrollment)')
    ylim([0.5 1])
    grid on
    ytickformat('%.2f')
    box on;
    hold off
    subplot(2,2,3)
    hold on
    set(gca,'TickLabelInterpreter','latex') 
    title('Female, Private')
    plot(xx,yy31,'Color', c(80,:),'LineWidth',1.2); 
    plot(xx,yy32,'Color', c(160,:),'LineWidth',1.2); 
    plot(xx,yy33,'Color', c(240,:),'LineWidth',1.2); 
    ylabel('Pr(enrollment)')
    ylim([0.5 1])
    grid on
    ytickformat('%.2f')
    box on;
    hold off
    subplot(2,2,4)
    hold on
    set(gca,'TickLabelInterpreter','latex') 
    title('Female, Public')
    plot(xx,yy41,'Color', c(80,:),'LineWidth',1.2); 
    plot(xx,yy42,'Color', c(160,:),'LineWidth',1.2); 
    plot(xx,yy43,'Color', c(240,:),'LineWidth',1.2); 
    ylabel('Pr(enrollment)')
    ylim([0.5 1])
    ytickformat('%.2f')
    box on;
    grid on
    legend({'2012','2011','2010'},'Location','southeast','Orientation','horizontal','Interpreter','latex')
    hL = legend({'2012','2011','2010'},'Location','southeast','Orientation','horizontal','Interpreter','latex');
    newPosition = [0.5 0 0.05 0.05];
    newUnits = 'normalized';
    set(hL,'Position', newPosition,'Units', newUnits);
    hold off
    % G25 Enrollment Probability By Gender/SES
    saveas(gcf, fullfile(pathOutput, 'figures', 'enrollmentandgrad_score_type.png'));
    exportgraphics(gcf, fullfile(pathOutput, 'figures','enrollmentandgrad_score_type.eps'),'BackgroundColor','none','ContentType','vector')

    
%% Figure 8: Model Fit - Selectivity
%(in Julia files)

%% Figure 9: BVP Impacts
%(in Julia files)
    
%% Figure 10: Impacts of Reducing Frictions (\alpha)
    
counterfactuals = readtable([pathData '/counterfactuals_big.csv']);
experiments=["scaleAlpha", "scaleAlpha2011"];
outcomes=["enroll", "welfare", "graduate"];
scales=["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"];
types=["Private Male", "Public Male", "Private Female", "Public Female"];
legends=["All On Platform","Exclude G8"];

%convert enrollment to indicator
for exp1=1:length(experiments)
for sca1=1:length(scales)
column1=experiments(exp1)+"_"+scales(sca1)+"_enroll"; 
counterfactuals.(column1)=counterfactuals.(column1)>0;
end
end

outcomes20112012=zeros(length(scales),length(outcomes),length(experiments));
for exp1=1:length(experiments)
for out1=1:length(outcomes)
for sca1=1:length(scales)
column1=experiments(exp1)+"_"+scales(sca1)+"_"+outcomes(out1);
outcomes20112012(sca1,out1,exp1)=mean(counterfactuals.(column1));
end
end
end
outcomes2011=outcomes20112012(:,:,2); outcomes2011(:,3)=outcomes2011(:,3)./outcomes2011(:,1);
outcomes2012=outcomes20112012(:,:,1); outcomes2012(:,3)=outcomes2012(:,3)./outcomes2012(:,1);


figure(111)
plot(0.9:-0.1:0,outcomes2012(:,1),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,1),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel("Enrollment Rate")
grid on
box on
hold off
set(gca,'TickLabelInterpreter','latex') 
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Enrollment.png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures','Counterfactual_ScaleAlpha_Enrollment.eps'),'BackgroundColor','none','ContentType','vector')

figure(112)
plot(0.9:-0.1:0,outcomes2012(:,2),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,2),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel("Welfare (1m Chilean Pesos)")
grid on
box on
hold off
set(gca,'TickLabelInterpreter','latex')
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Welfare.png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures','Counterfactual_ScaleAlpha_Welfare.eps'),'BackgroundColor','none','ContentType','vector')

figure(113)
plot(0.9:-0.1:0,outcomes2012(:,3),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,3),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel("Graduation Rate")
grid on
box on
ylim([0.5 0.535])
hold off
set(gca,'TickLabelInterpreter','latex')
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Graduation.png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures','Counterfactual_ScaleAlpha_Graduation.eps'),'BackgroundColor','none','ContentType','vector')

%% Figure 11: Welfare Impacts of Reducing Frictions
for type=1:4

outcomes20112012=zeros(length(scales),length(outcomes),length(experiments));
for exp1=1:length(experiments)
for out1=1:length(outcomes)
for sca1=1:length(scales)
column1=experiments(exp1)+"_"+scales(sca1)+"_"+outcomes(out1);
outcomes20112012(sca1,out1,exp1)=mean(counterfactuals.(column1)(counterfactuals.type==type));
end
end
end
outcomes2011=outcomes20112012(:,:,2); outcomes2011(:,3)=outcomes2011(:,3)./outcomes2011(:,1);
outcomes2012=outcomes20112012(:,:,1); outcomes2012(:,3)=outcomes2012(:,3)./outcomes2012(:,1);


figure(201)
plot(0.9:-0.1:0,outcomes2012(:,1),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,1),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel(['Enrollment Rate | ' + types(type)])
grid on
box on
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
legend({legends(1),legends(2)},'Location','southeast','Interpreter','latex')
legend('boxoff') 
xlim([0 0.9])
if type == 3
    ylim([0.66 0.685])
else
end
hold off
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Enrollment" + "_type" +num2str(type) + ".png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Enrollment" + "_type" +num2str(type) + ".eps"),'BackgroundColor','none','ContentType','vector')

figure(202) %Figure 11 contains "figure(202)" for all 4 types
plot(0.9:-0.1:0,outcomes2012(:,2),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,2),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel(['Welfare | ' + types(type)])
grid on
box on
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
legend({legends(1),legends(2)},'Location','southeast','Interpreter','latex')
legend('boxoff')
xlim([0 0.9])
hold off
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Welfare" + "_type" +num2str(type) + ".png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Welfare" + "_type" +num2str(type) + ".eps"),'BackgroundColor','none','ContentType','vector')

figure(203)
plot(0.9:-0.1:0,outcomes2012(:,3),'Color', c(80,:),'LineWidth',3);
hold on
plot(0.9:-0.1:0,outcomes2011(:,3),'Color', c(160,:),'LineWidth',3);
xlabel("Fraction Reduction in Frictions")
ylabel(['Graduation Rate | ' + types(type)])
grid on
box on
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
legend({legends(1),legends(2)},'Location','southeast','Interpreter','latex')
legend('boxoff')
xlim([0 0.9])
if type == 1
    ylim([0.51 0.535])
elseif type == 2
    ylim([0.39 0.43])
elseif type == 3
	ylim([0.675 0.7])
elseif type == 4
	ylim([0.56 0.6])
else
end
hold off
saveas(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Graduation" + "_type" +num2str(type) + ".png"));
exportgraphics(gcf, fullfile(pathOutput, 'figures',"Counterfactual_ScaleAlpha_Graduation" + "_type" +num2str(type) + ".eps"),'BackgroundColor','none','ContentType','vector')
end


%% Figure 12: Utility loss of removing options ordered by selectivity

counterfactualstype = readtable([pathData '/counterfactuals_bytype.csv']);
keepcolumns= {'inds','benchmark_welfare','dropSelectivityBin1_welfare','dropSelectivityBin2_welfare','dropSelectivityBin3_welfare','dropSelectivityBin4_welfare',...
    'dropSelectivityBin5_welfare','dropSelectivityBin6_welfare','dropSelectivityBin7_welfare','dropSelectivityBin8_welfare', ...
    'dropSelectivityBin9_welfare','dropSelectivityBin10_welfare'};
dropSelectivity = counterfactualstype(:,keepcolumns);
types=["Private Male", "Public Male", "Private Female", "Public Female"];

% All
rows = [1:5:130];
dropSelectivity = dropSelectivity(rows,:);
avgs = mean(table2array(dropSelectivity(:,[2:12])));
% Relative to benchmark
avgs(2:11) = avgs(2:11)-avgs(1);


figure(13)
h = bar([1:10],avgs(2:11),'FaceColor', c(40,:),'EdgeColor', c(20,:),'FaceAlpha', 0.8);
ylabel('Welfare Loss')
xlabel('Decile of selectivity')
grid on
box on
xlim([0 12])
yline(avgs(3),'--')
yline(avgs(11),'--')

r = -(avgs(3)-avgs(11))/avgs(3);

% These define the placement and size of the brace
y1 = avgs(3);
y2 = avgs(11);
x = 10.8;
width = 0.2;
% Make some x-data and y-data
line_y = y1  + [0, 0.02, 0.48, 0.5, 0.52, 0.98, 1]*(y2-y1);
line_x = x  + [0, 0.25, 0.25, 0.6, 0.25, 0.25, 0]*width;
% Draw the brace and some text
line(line_x, line_y, 'Color', 'k')
text(11,-0.17,[num2str(round(r*100)) '%']);
text(11,-0.18,'of D2');
ytickformat('%.2f')
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
hold off
saveas(gcf, fullfile(pathOutput, 'figures','utility_loss_selectivity.png'));
exportgraphics(gcf, fullfile(pathOutput, 'figures','utility_loss_selectivity.eps'),'BackgroundColor','none','ContentType','vector')

%% Figure 13: Change in Welfare, Enrollment and Graduation Rates by Type
counterfactualstype.Properties.VariableNames
keepcolumns= {'inds','benchmark_welfare','benchmark_graduation','benchmark_enroll','asif2011_welfare','asif2011_graduation','asif2011_enroll'};

changes = counterfactualstype(:,keepcolumns);
types=["Male Private", "Male Public", "Female Private", "Female Public"];

% All
rowsMPriv = [2:5:130];
rowsMPubl = [3:5:130];
rowsFPriv = [4:5:130];
rowsFPubl = [5:5:130];

MPriv = changes(rowsMPriv,:);
MPubl = changes(rowsMPubl,:);
FPriv = changes(rowsFPriv,:);
FPubl = changes(rowsFPubl,:);


MPriv = mean(table2array(MPriv(:,[2:7])));
MPriv = [MPriv(1) - MPriv(4) MPriv(2) - MPriv(5) MPriv(3) - MPriv(6)];

MPubl = mean(table2array(MPubl(:,[2:7])));
MPubl = [MPubl(1) - MPubl(4) MPubl(2) - MPubl(5) MPubl(3) - MPubl(6)];

FPriv = mean(table2array(FPriv(:,[2:7])));
FPriv = [FPriv(1) - FPriv(4) FPriv(2) - FPriv(5) FPriv(3) - FPriv(6)];

FPubl = mean(table2array(FPubl(:,[2:7])));
FPubl = [FPubl(1) - FPubl(4) FPubl(2) - FPubl(5) FPubl(3) - FPubl(6)];

changeswelfare = [MPriv(1) MPubl(1) FPriv(1) FPubl(1)];
changesgraduation = 100*[MPriv(2) MPubl(2) FPriv(2) FPubl(2)];
changesenroll = 100*[MPriv(3) MPubl(3) FPriv(3) FPubl(3)];

figure(141)
hold on
b = bar(1,changeswelfare);
grid on
set(b(1),'FaceColor',c(10,:),'EdgeColor',c(10,:))
set(b(2),'FaceColor',c(70,:),'EdgeColor',c(70,:))
set(b(3),'FaceColor',c(110,:),'EdgeColor',c(110,:))
set(b(4),'FaceColor',c(150,:),'EdgeColor',c(150,:))
legend(types,'Location','northwest','Interpreter','latex')
box on
ylabel('Expected change in welfare')
xlabel('As if 2012 vs. As if 2011')
ytickformat('%.2f')
set(gca,'FontSize', 14)
xticks({})
set(gca,'TickLabelInterpreter','latex')
hold off
saveas(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_welfare.png'));
laprint(141,fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_welfare'))
exportgraphics(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_welfare.eps'),'BackgroundColor','none','ContentType','vector')

figure(142)
hold on
b = bar(1,changesgraduation);
grid on
set(b(1),'FaceColor',c(10,:),'EdgeColor',c(10,:))
set(b(2),'FaceColor',c(70,:),'EdgeColor',c(70,:))
set(b(3),'FaceColor',c(110,:),'EdgeColor',c(110,:))
set(b(4),'FaceColor',c(150,:),'EdgeColor',c(150,:))
legend(types,'Location','northwest','Interpreter','latex')
box on
ylabel('Expected change in graduate (\%)')
xlabel('As if 2012 vs. As if 2011')
xticks({})
ytickformat('%.1f')
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
hold off
saveas(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_graduate.png'));
exportgraphics(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_graduate.eps'),'BackgroundColor','none','ContentType','vector')


figure(143)
hold on
b = bar(1,changesenroll);
grid on
set(b(1),'FaceColor',c(10,:),'EdgeColor',c(10,:))
set(b(2),'FaceColor',c(70,:),'EdgeColor',c(70,:))
set(b(3),'FaceColor',c(110,:),'EdgeColor',c(110,:))
set(b(4),'FaceColor',c(150,:),'EdgeColor',c(150,:))
legend(types,'Location','northwest','Interpreter','latex')
box on
ylabel('Expected change in enrollment (\%)')
xlabel('As if 2012 vs. As if 2011')
xticks({})
ytickformat('%.0f')
ax = gca;
ax.YRuler.Exponent = 0;
set(gca,'FontSize', 14)
set(gca,'TickLabelInterpreter','latex')
hold off
saveas(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_enroll.png'));
exportgraphics(gcf, fullfile(pathOutput, 'figures','Bars_benchmark_asif2011_enroll.eps'),'BackgroundColor','none','ContentType','vector')
